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Chapter 5 

Finite Element Method 



Guillaume Demesy, Frederic Zolla, Andre Nicolet, and Benjamin Vial 

Aix-Marseille Universite, Ecole Centrale Marseille, Institut Fresnel, 
13397 Marseille Cedex 20, France 

iguillaum e.demesy @fresnel.fr^ 

5.1 Introduction 

Finite element methods (FEM) represent a very general set of techniques to approximate solu- 
tions of partial derivative equations. Their main advantage lies in their ability to handle arbitrary 
geometries via unstructured meshes of the domain of interest: The discretization of oblic ge- 
ometry edges is natively built in. Finite Element Methods have been widely developed in many 
areas of physics and engineering: mechanics, thermodynamics. . . 

But until the early 80's, two major drawbacks prevented them from being used in electro- 
magnetic problems. On the one hand, existing nodal element basis did not satisfy the physical 
(dis)continuity of the vector fields components and lead to spurious solutions [IJ. On the other 
hand, there was no proper way to truncate unbounded regions in open wave problems. 

These two major limitations were both overcome in the early 80's: Vector elements have 
been developed by Nedelec 0121, and Perfectly Matched Layers (PMLs) were discovered by 
Berenger [4J. Since then, it has been shown that PMLs could be described in the general frame- 
work of transformation optics [|5l[6l|7l|8]|. 

All the mathematical and computational ingredients now exist and the goal of this chapter 
is to show how to combine them to implement a general 3D numerical scheme adapted to 
gratings using Finite Elements. In fact, we are now facing the physical difficulties inherent to 
the infinite spatial characteristics of the grating problem, whereas the computation domain has 
to be bounded in practice: (i) Both the superstrate and the substrate are infinite regions, (ii) 
there is an infinite number of periods and, last but not least, (iii) the sources of the incident field 
(a plane wave) are located in the superstrate at an infinite distance from the grating. 

In this chapter, the infinite extension of the superstrate and substrate is addressed using 
cartesian PMLs. In the framework of transformation optics, we demonstrate that Berenger's 
original PMLs can be extended to the challenging numerical cases of grazing incidence in order 
to deal with extreme oblic incidences or configurations near Wood's anomalies. The second is- 
sue of infinite number of period can be addressed via Bloch conditions. Finally, we are dealing 
with the distant plane wave sources through an equivalence of the diffraction problem with a ra- 
diation one whose sources are localized inside the diffractive element itself. The unknown field 
to be approximated using Finite Elements is a radiated field with sources inside the computation 
box and allows to retrieve easily the total field with the plane wave source. 

In a first section, we derive and implement this approach in the so-called 2D non-conical, 
or scalar, case. We are dealing with the infinite issues rigorously in both TE and TM polarization 
cases. It results in a radiation problem with sources locahzed in the diffractive element itself. 
We mathematically split the whole problem into two parts. The first one consists in the classical 
calculation of the total field solution of a simple interface. The second one amounts to looking 
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for a radiated field with sources confined within the diffractive obstacles and deduced from 
the first elementary problem. From this viewpoint, the later radiated field can be interpreted 
as an exact perturbation of the total field. We show that our approach allows to tackle some 
kind of anisotropy without increasing the computational time or resource. Through a battery of 
examples, we illustrate its independence towards the geometry of the diffractive pattern. Finally, 
we present an Adaptative PML able to tackle grazing incidences or configurations near Wood's 
anomaly. 

In a second section, we extend this approach to the most general configuration of vector 
diffraction by crossed gratings embedded in arbitrary multilayered stack. The main advantage 
of this method is, again, its complete independence towards the shape of the diffractive element, 
whereas other methods often require heavy adjustments depending on whether the geometry of 
the groove region presents oblique edges. This approach combined with the use of second order 
edge elements allows us to retrieve the few numerical academic examples found in the literature 
with an excellent accuracy. Furthermore, we provide a new reference case combining major 
difficulties: A non trivial toroidal geometry together with strong losses and a high permittivity 
contrast. Finally, we discuss computation time and convergence as a function of the mesh 
refinement as well as the choice of the direct solver. 



5.2 Scalar diffraction by arbitrary mono-dimensional gratings : a Finite Element for- 
mulation 

5.2.7 Set up of the problem and notations 

We denote by x, y and z, the unit vectors of the axes of an orthogonal coordinate system Oxyz. 
We deal only with time-harmonic fields; consequently, the electric and magnetic fields are rep- 
resented by the complex vector fields E and H, with a time dependance in exp(— / 0)^). 

Besides, in this chapter, we assume that the tensor fields of relative permittivity e and 
relative permeability ji can be written as follows: 

/ l^xx \ 

and Ai= Ma liyy \ , (5.1) 

where E^x.^a," l^zz possibly complex valued functions of the two variables x and y and 
where £a (resp. jla) represents the conjugate complex of £a (resp. jia). These kinds of materials 
are said to be z-anisotropic. It is of importance to note that with such tensor fields, lossy 
materials can be studied (the lossless materials correspond to tensors with real diagonal terms 
represented by Hermitian matrices) and that the problem is invariant along the z-axis but the 
tensor fields can vary continuously (gradient index gratings) or discontinuously (step index 
gratings). Moreover we define ko := co/c. 



^xx 









£yy 











£zz 



The gratings that we are dealing with are made of three regions (See Fig. 5.1 ). 



The superstratum {y > h^) which is supposed to be homogeneous, isotropic and lossless 
and characterized solely by its relative permittivity and its relative permeability fi^ 
and we denote k^ ko y^£+/i+ 

The substratum (y < 0) which is supposed to be homogeneous and isotropic and there- 
fore characterized by its relative permittivity £~ and its relative permeability /i~ and we 
denote k~ := ko ^j£~\i~ 
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The groove region (0 <y < h^) which can be heterogeneous and z-anisotropic and thus 
characterized by the two tensor fields £^{x^y) and ii^{x^y). It is worth noting that the 
method does work irrespective of whether the tensor fields are piecewise constant. The 
groove periodicity along x-axis will be denoted d. 
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Fig. 5.1: Sketch and notations of the grating studied in this section. 



This grating is illuminated by an incident plane wave of wave vector k+ = ax — j8+y = 
(sin^ox — cos^oy). whose electric field (TM case) ( resp. magnetic field (TE case)) is lin- 
early polarized along the z-axis: 

eO = AO exp(/k; • r) z (resp. RO = exp(/k; • r) z) , (5.2) 

where (resp. A^) is an arbitrary complex number and r — {x^yY . In this section, a plane 
wave is characterized by its wave- vector denoted kj^'^^^ The subscript p (resp. c) stands 
for "propagative" (resp. "counter-propagative"). The superscript + (resp. — ) refers to the 
associated wavenumber (resp. k~), and indicates that we are dealing with a plane wave 
propagating in the superstrate (resp. substrate). The magnetic (resp. electric) field derived from 
(resp. H^) is denoted (resp. E|J^) and the electromagnetic field associated with the 
incident field is therefore denoted (E^,H^) which is equal to (E^,H^) (resp. (E^,H^)). 

The diffraction problem that we address consists in finding Maxwell equation solutions 
in harmonic regime i.e. the unique solution (E,H) of: 

curl E = / 0) ^0 M H (5.3a) 
curlH = -i(0£{)£E (5.3b) 

such that the diffracted field (E^, H^) := (E - E^, H - H^) satisfies an Outgoing Waves Condi- 
tion (O.W.C. [9]) and where E and H are quasi-periodic functions with respect to the x coordi- 
nate. 
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5.2.2 Theoretical developments of the method 

5.2.2.1 Decoupling of fields and z-anisotropy 

We assume that 5_{x,y) is a z-anisotropic tensor field (5;^^ = = 5^ = S^y — 0). Moreover, the 
left upper matrix extracted from 5 is denoted 5, namely: 




(5.4) 



For z-anisotropic materials, in a non-conical case, the problem of diffraction can be split into 
two fundamental cases (TE case and TM case). This property results from the following equality 
which can be easily derived: 

- curl {dr^ curl (^^z)) = div (|^/det(|)Vi^) z , (5.5) 

where is a function which does not depend on the z variable. Relying on the previous equality, 
it appears that the problem of diffraction in a non conical mounting amounts to looking for 
an electric (resp. magnetic) field which is polarized along the z-axis ; E = ^(x,_y)z (resp. 
H = h{x^y) z). The functions e and h are therefore solutions of similar differential equations: 

^^^^{u) := div (l^Vw) +klxu = (5.6) 

with 

u = e, £ = £^/det(£), X = ^zz ^ (5.7) 

in the TM case and 

u = h, £ = |^/det(|), Z = Mzz, (5.8) 

in the TE case. 



5.2.2.2 Boiling down the diffraction problem to a radiation one 



In its initial form, the diffraction problem summed up by Eq. ( |5.6| ) is not well suited to the 
Finite Element Method. In order to overcome this difficulty, we propose to split the unknown 
function u into a sum of two functions u\ and the first term being known as a closed form 
and the latter being a solution of a problem of radiation whose sources are localized within the 
obstacles. 



We have assumed that outside the groove region (cf. Fig. 5.1), the tensor field ^ and 

the function x are constant and equal respectively to ^~ and ;f ~ in the substratum (y < 0) and 

equal respectively to and x^ in the superstratum (y > h^). Besides, for the sake of clarity, 
the superstratum is supposed to be made of an isotropic and lossless material and is therefore 
solely defined by its relative permittivity and its relative permeability which leads to: 

1^+ = ^ Id2 and x^ = in TE case (5.9) 

or 

1^+ = — Id2 and x^ = in TM case, (5.10) 
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where Id2 is the 2 x 2 identity matrix. With such notations, ^ and x ^re therefore defined as 
follows: 

r y>^' ( for y> ¥ 

^{x,y):={ l'{x,y) for ¥>y>0 ^ ^^^^y) ) ^g^^^y) for ¥>y>0 (5.11) 
r for J < [ X~ for J < . 

It is now apropos to introduce an auxiliary tensor field ^ ^ and an auxiliary function Xi • 

{^^ for y>0 ( y+ for ^>o 

1- for ,<0 •^'(^■^'-{z- /<0 <='^> 

these quantities corresponding, of course, to a simple plane interface. Besides, we introduce the 
constant tensor field which is equal to everywhere and a constant scalar field Xo which is 

equal to x^ everywhere. Finally, we denote uq the function which equals the incident field u^^^ 
in the superstratum and vanishes elsewhere (see Fig. |5.1[ ): 

. . r ^i^^ for y>¥ 
-o(^'^)^=| for y<h^ ^^-13^ 

We are now in a position to define more precisely the diffraction problem that we are 
dealing with. The function u is the unique solution of: 

J5f^^(i^) = 0, such that u'^ := u — uq satisfies an O.W.C. (5.14) 

In order to reduce this diffraction problem to a radiation problem, an intermediate function is 
necessary. This function, called wi, is defined as the unique solution of the equation: 

M (^i) ~ ^ ' ^^^^ ^^^^ := u\—uo satisfies an O.W.C. (5.15) 

The function ui corresponds thus to an annex problem associated to a simple interface and can 
be solved in closed form mdfrom now on is considered as a known function. As written above, 
we need the function ^2 which is simply defined as the difference between u and u\: 

U2'=u — ui=u^ — Ui. (5.16) 

The presence of the superscript d is, of course, not irrelevant: As the difference of two diffracted 
fields, the O.W.C. of ^2 is guaranteed (which is of prime importance when dealing with PML 
cf. |5.2.2.4| ). As a result, the Eq. ( |5.14| ) becomes: 

^i,z(^2) = -^£,z(^i)' (5.17) 

where the right hand member is a scalar function which may be interpreted as a known source 
term —S^\{x^y) and the support of this source is localized only within the groove region. To 
prove it, all we have to do is to use Eq. ( |5.15| ): 
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Now, let us point out that the tensor fields ^ and ^ ^ are identical outside the groove region and 

the same holds for x arid Xi • The support of S^i is thus localized within the groove region as 
expected. It remains to compute more explicitly the source term S^i . Making use of the linearity 
of the operator J5f and the equality ui = u^ + uoA^^^ source term can be split into two terms 



where 



and 



(5.19) 
(5.20) 
(5.21) 



Now, bearing in mind that uq is nothing but a plane wave uq = exp(/k+ • r) (with k+ = ax - 
j8+y), it is sufficient to give Vwq = /k+ uq for the weak formulation associated with Eq. (5.17): 



^0 = {/div -i) exp(/k; -r)] +kl (x^-x) exp(/k; -r) 



(5.22) 



The same holds for the term associated with the diffracted field. Since, in the superstrate, we 
have of course u'f = p exp(«k+ • r) with k+ = ax + /3+y, 

^/ = p{/div[(£+-£)k+exp(/k+T)] +kl{x+-x)^MiK-r)} , (5-23) 
where p is simply the complex reflection coefficient associated with the simple interface: 

jS^intheTM case 

in the TE case 



P = 



1+ . 



with = 



(5.24) 



5.2.2.3 Quasi-periodicity and weak formulation 

The weak formulation follows the classical lines and is based on the construction of a weighted 



residual of Eq. ( |5.6| ), which is multiplied by the complex conjugate of a weight function and 
integrated by part to obtain: 



^^^^{u,u) = - j y^VuyVu' + klxuu' + y^VuyndS (5.25) 

The solution u of the weak formulation can therefore be defined as the element of the space 
L^(curl,(i,a) of quasiperiodic functions (i.e. such that u{x^y) — u#{x^y)e^^^ with u#{x^y) = 
u#{x + d,y), a J-periodic function) of L^(curl) on such that: 

^^^^{u,u')=0 Vi/' gL^ (curl, J, a). (5.26) 

As for the boundary term introduced by the integration by part, it can be classically set to zero 
by imposing Dirichlet conditions on a part of the boundary (the value of u is imposed and the 
weight function can be chosen equal to zero on this part of the boundary) or by imposing 
homogeneous Neumann conditions ((^ Vw) • n = on another part of the boundary (and u is 
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Fig. 5.2: Quasi-periodicity of the field and sample of a d-periodic mesh. 



therefore an unknown to be determined on the boundary). A third possibihty are the so-called 
quasi-periodicity conditions of particular importance in the modeling of gratings. 

Denote by F/ and Tr the lines parallel to the };-axis delimiting a cell of the grating (see 
Fig. |5.2| ) respectively from its left and right neighbor cell. Considering that both u and are in 
L^(curl, d, a), the boundary term for F/ U F^ is 



/ Vu]-ndS= [ 4^-^"^^ V(w#^+^'^^)') . n d5 = 

/ ^# ( <^ ( Vw# + iau#x) ] • n dS = , 



/r/ur, 

because the integrand u'^ (^^ {Vu# + iau#x)^ • n is periodic along x and the normal n has opposite 

directions on F/ and F^ so that the contributions of these two boundaries have the same absolute 
value with opposite signs. The contribution of the boundary terms vanishes therefore naturally 
in the case of quasi-periodicity. 

The finite element method is based on this weak formulation and both the solution and the 
weight functions are classically chosen in a discrete space made of linear or quadratic Lagrange 
elements, i.e. piecewise first or second order two variable polynomial interpolation built on a 



triangular mesh of the domain (cf. Fig. |5.3a| ). Dirichlet and Neumann conditions may be used 
to truncate the PML domain in a region where the field (transformed by the PML) is negligible. 
The quasi-periodic boundary conditions are imposed by considering the u as unknown on F/ 
(in a way similar to the homogeneous Neumann condition case) while, on F^, u is forced equal 
to the value of the corresponding point on F/ (i.e. shifted by a quantity —d along x) up to the 
factor e^^^. The practical implementation in the finite element method is described in details in 
MM 



5.2.2.4 Perfectly Matched Layer for z-anisotropic materials 



The main drawback encountered in electromagnetism when tackling theory of gratings through 
the finite element method is the non-decreasing behaviour of the propagating modes in super- 
stratum and substratum (if they are made of lossless materials): The PML has been introduced 
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by [m in order to get round this obstacle. The computation of PML designed for z-anisotropic 
gratings is the topic of what follows. 

In the framework of transformation optics, a PML may be seen as a change of coordinate 
corresponding to a complex stretch of the coordinate corresponding to the direction along which 
the field must decay [[III [III (Ml. Transformation optics have recently unified various techniques 
in computational electromagnetics such as the treatment of open problems, helicoidal geome- 
tries or the design of invisibility cloaks ([T5]). These apparently different problems share the 
same concept of geometrical transformation, leading to equivalent material properties. A very 
simple and practical rule can be set up ([lOJ): when changing the coordinate system, all you 
have to do is to replace the initial materials properties £ and ji by equivalent material properties 
£y and lis given by the following rule: ~ 

= j-i^j-Tdet(J) and = J" V J"^det(J), (5.27) 

where J is the Jacobian matrix of the coordinate transformation consisting of the partial deriva- 
tives of the new coordinates with respect to the original ones (J~^ is the transposed of its 
inverse). 

In this framework, the most natural way to define PMLs is to consider them as maps on a 
complex space C^, which coordinate change leads to equivalent permittivity and permeability 
tensors. We detail here the different coordinates used in this section. 

• {x,y,z) are the cartesian original coordinates. 

• {xs,ys,Zs) are the complex stretched coordinates. A suitable subspace F c is chosen 
(with three real dimensions) such that {xs^ys^Zs) are the complex valued coordinates of a 
point on T (e.g. x = 9^e(x^), y = yit{ys). z = 9^e(z^)). 

• {xc^yc^Zc) are three real coordinates corresponding to a real valued parametrization of 

rcc^. 

We use rectangular PMLs ([fT2ll) absorbing in the _y-direction and we choose a diagonal 
matrix J = 6idig(\ ^ Sy{y) ^ 1), where Sy{y) is a complex- valued function of the real variable y, 
defined by: 

ys{y)= Fsy{y')dy', (5.28) 
The expression of the equivalent permittivity and permeability tensors are thus: 

/ SyEjcx £a \ / Syll^^ JTa \ 

^= Sy^Eyy and ^^5= lla Sy^llyy . (5.29) 

\ Sy£^z I ~ \ Syll^z I 



Note that the equivalent medium has the same impedance than the original one as £ an 
are transformed in the same way, which guarantees that the PML is perfectly reflectionless. ~ 
Now, let us define the so-called substituted field Fs — (E^,H^), solution of Eqs. (5.3) with 
t,—t,s and X — Xs' It turns out that Fs equals the field F in the region y^ <y <y^ (with y^ — —h~ 
and y = h^ + h^ , see Fig. 5.3a), provided that Sy{y) = 1 in this region. The main feature of this 
latest field Fs is the remarkable correspondence with the first field F ; whatever the function 
Sy provided that it equals 1 for < _y < _y^, the two fields F and Fs are identical in the region 
y^ <y < j^[8J. In other words, the PML is completely reflection-less. In addition, for complex 
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valued functions Sy {Zm{sy} strictly positive in PML), the field Fs converges exponentially 



towards zero (as y tends to ±00, cf. Fig. |5.3c| and |5.3d| ) although its physical counterpart F 
does not. Note that in Fig. 



5.3d 



the value of the computed radiated field ^2 ^^^h extreme 
boundary of the PMLs is at least 10~^ weaker than in the region of interest. As a consequence, 
Fs is of finite energy and for this substituted field a weak formulation can be easily derived 
which is essential when dealing with Finite Element Method. 

Still remains to give a suitable function Sy. Let us consider the complex coordinate map- 
ping y{yc), which is simply defined as the derivative of the stretching coefficient Sy{y) with 
respect to yc. With simple stretching functions, we can obtain a reliable criterion upon proper 
fields decay. A classical choice is: 



if y<y^ 

if y^ <y<y' 

ify>f 



(5.30) 



where = ^ + are complex constants with ^' '^ > 0. 

In that case, the complex valued function y{yc) defined by Eq. (5.28) is explicitly given by: 



y{yc) 



'/ + C-(}^c-/) ifjc</ 

yc if y^ <yc<y^ ^ 

y' + ^^{yc-y') '^^yc>y' 



(5.31) 



Finally, let us consider a propagating plane wave in the substratum Un{x^y) := exp(/(ax — 
j8~_y)). Its expression can be rewritten as a function of the stretched coordinates in the PML as 
follows: 

u'^{xc,yc) := Un{x{xc)Myc)) = e'^-^e-'^n{^^^-{yc-n (5.32) 

The behavior of this latest function along the yc direction is governed by the function lJ^^{y^ \— 
e-iPnC-yc^ Letting i8^'- :=^e{p-}, pH^- :=3m{i8-}, C^- :=MC"} and C'^' :=3m{C-}, 
the non-oscillating part of the function U^^{yc) is given by exp (^{Pn~ + C,'~)y^ . 
Keeping in mind that and/or are positive numbers, the function decreases expo- 



nentially towards zero as yc tends to — oo (Fig. 5.3d) provided that ^~ belongs to C+ := {z G 
C,9t^{z} > 0, and 3m{z} > 0}. In the same way, it can be shown that ^+ belongs to C+. 
Let us conclude this section with two important remarks: 

1. Practical choice of PML parameters. As for the complex stretch parameters, setting 
= 1 + / is usually a safe choice. For computational needs, the PML has t o be t runcated 



and the other constitutive parameter of the PML is its thickness h (see Fig. |5.3a| ). Setting 
/i^ = Ao/ ^/e^ leads to a PML thick enough to "absorb" all incident radiation. These 
specific values will be used in the sequel, unless otherwise specified. 

2. Special cases. The reader will notice that a configuration where is a very weak 
positive number compared to with (this is precisely the case of a plane wave at 
grazing incidence on the bottom PML) leads to a very slow exponential decay of W^. 

In such a case, close to so-called Wood's anomalies or at extreme grazing incidences. 



classical PML fail. We will address this tricky situation extensively in Section 5.2.4 
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(a) Computational 
domain and its five 
constituent regions. 





(b) Coarse triangle (c) Radiated field: (d) Radiated field: 

meshing of the 3ie{u2}inV/m log(|w2l) 

cell Q,. Maximum 

element side size: 

X/{2Vi) 



Fig. 5.3: Example of computation of the radiated field U2 (TM case). 



5.2.2.5 Synthesis of the method 

In order to give a general view of the method, all information is collected here that is neces- 
sary to set up the practical Finite Element Model. First of all, the computation domain (cf. 
Fig. |5.3a| ) corresponds to a truncated cell of the grating which is a finite rectangle divided into 



five horizontal layers. These layers are respectively from top to bottom upper PML, the su- 
perstratum, the groove region, t he sub stratum, and the lower PML. The unknown field is the 
scalar function uf defined in Eq. ( 5.16[ ). Its finite element approximation is based on the second 



Lagrange elements built on a triangle meshing of the cell (cf. Fig. |5.3b| ). A complex algebraic 
system of linear equations is constructed via the Galerkin weighted residual method, i.e. the set 
of weight functions is chosen as the set of shape functions of interpolation on the mesh MOil . 



In region 1 (upper PML, see Fig. |5.3a| ), 

^^^^^^{uiu')=0, (5.33) 

=s 

with ^ and depending on the equivalent anisotropic properties of the PML given by 
Eq. ^J}, Eq. (Q and Eqs. ^29^ . 

In region 2 (superstratum), 

^^^^^4uiu')=0, (5.34) 

with and depending on the homogeneous isotropic properties of the superstratum 
given'by Eq. ([5J]), Eq. (|5]8]), Eq. and Eq. ( [5?T0| ). 
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In region 3 (groove region), 



(5.35) 



with and depending on the heterogeneous possibly anisotropic properties given by 
Eq. ^J}, Eq. Eq. ^AA\ and given by Eq. ( |5l9l ) , Eq. ( |5^ , Eq. ( |5^ and 
Eq. (|5:24l). 



In region 4 (substratum). 



?^-.^_(«f,«')=0, 



(5.36) 



with £, and ;^ depending on the homogeneous isotropic properties of the substratum 
givenT)y Eq. ^J), Eq. (|5j8]), Eq. (|5]9]) and Eq. ( |5l0l ). 



In region 5 (lower PML), 



(5.37) 



with ^ and depending on the equivalent anisotropic properties of the PML given by 
Eq. ^J}, Eq. and Eqs. ^529^ . 



5.2.2.6 Energy balance: Diffraction efficiencies and absorption 



The rough result of the FEM calculation is the complex radiated field i^2- Usi ng Eg . ( 5.16| ), 
it is straightforward to obtain the complex diffracted field solution of Eq. ( 5.14 ) at each 
point of the bounded domain. We deduce from the diffraction efficiencies with the following 
method. The superscripts + (resp. ~) correspond to quantities defined in the superstratum (resp. 
substratum) as previously. 

On the one hand, since is quasi-periodic along the x-axis , it can be expanded as a 
Rayleigh expansion (see for instance 



where 



4{y) = 



1 rd/2 



2n 

u" (x, y)e '""^ dx with a„ = a + — n 



(5.38) 



(5.39) 



d J-d/2 d 
On the other hand, introducing Eq. ( |5.38[ ) into Eq. ( |5.6[ ) leads to the Rayleigh coefficients: 

i+(y) = rne'Pny ^ane-'Pny for y>hs 



<{y) 



u-{y)=tne-'^-y + bne'^-y for y<Q 



withj8„^ 



(5.40) 



For a temporal dependance in e the O.W.C. imposes an = bn = 0. Combining Eq. ( |5.39| ) 
and ( |5.40| ) at a fixed yo altitude leads to: 



r„ = i/" M'^(x,yo)e-'^""'^+'^«^^o)dx for yo>hS 

J-d/2 
rd/2 

tn = \ u^x^yo) dx for < 



(5.41) 
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We extract these two coefficients by trapezoidal numerical integration along x from a cutting of 
the previously calculated field map at _yo- It is well known that the mere trapezoidal integration 
method is very efficient for smooth and periodic functions (integration on one period) [16J. 
Now the restriction on a horizontal straight line crossing the whole cell in homogeneous media 
(substratum and superstratum) is of class. From a numerical point of view, it appears that the 
interpolated approximation of the unknown function, namely U2 preserves the good behaviour 
of the numerical computation of these integrals. From this we immediately deduce the reflected 
and transmitted diffracted efficiencies of propagative orders (Tn and Rn) defined by: 

«+ 

Rn'= rnrn^ foY yo > ¥ f 1 in the TM case 

with = I (5.42) 
T,:=tnrnff for 3;o<0 I in the TE case 

This calculation is performed at several different yo altitudes in the superstratum and the sub- 
stratum, and the mean value found for each propagative transmitted or reflected diffraction order 
is presented in the numerical experiments of the following section. 

Normalized losses Q can be obtained according to Poynting's theorem through the straight- 
forward computation of the following ratio: 

9?4EoxH^}.nd/ 



Q := , (5.43) 



L 



The numerator in Eq. ( |5.43| ) clarifies losses in Watts by period of the considered grating and are 
computed by integrating the Joule effect losses density over the surface S of the lossy element. 
The denominator normalizes these losses to the incident power, i.e. the time-averaged incident 
Poynting vector flux across one period (a straight line L of length d in the superstrate parallel to 
Ox, whose normal oriented along decreasing values of y is denoted n). 



Finally, combining Eqs. ( |5.42| ) and Eq. ( |5.43| ), a self consistency check of the whole 



numerical scheme consists in comparing the quantity B: 

B'.^Y.^n+Y^Rm + Q (5.44) 



to unity. In Eq. ( |5.44| ), the summation indexed by n (resp. m) corresponds to the sum over the 



efficiencies of all transmitted (resp. reflected) propagative diffraction orders in the substrate 
(resp. superstrate). We give interpretations and concrete examples of such numerical energy 
balances over non trivial grating profiles in sections [5.2.3.2| and |5.2.3.3 , 



5.2.3 Numerical experiments 

5.2.3.1 Numerical validation of the method 

We can refer to [17J in order to test the accuracy of our method. The studied grating is isotropic, 
since we lack numerical values in the literature in anisotropic cases. We compute the following 
problem (cf. Fig.|5.4|), as described in [18J and ifTTl . The wavelength of the plane wave is set to 



1 iim and is incoming with an angle of ;r/6 with respect to the normal to the grating. 

We present the Rq efficiency (cf. Table |5.1| ) in both cases of polarization versus the mesh 
refinement. So we have a good agreement to the reference values, and the accuracy reached is 
independent from the polarization case. 
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500 nm 

^ ^ superstrate 







1 \mi 


















1 |Lim 


Wate = - 44.9757 + 2.9524/ 



Fig. 5.4: Rectangular groove grating: This pattern is repeatedly set up with a period d = 1 jlm. This 
grating has been studied by U7J and is one of our points of reference 



Maximum element size 


r,TH 

^0 


r,TM 

^0 




0.7336765 


0.8532342 




0.7371302 


0.8456592 




0.7347466 


0.8482817 




0.7333739 


0.850071 




0.7346569 


0.8494844 


V(14V^) 


0.7341944 


0.8483238 




0.7342714 


0.8484774 


Result given by |17|| 


0.7342789 


0.8484781 



Tab. 5.1: Reflected efficiencies versus mesh refinement. Note that the efficiencies are properly com- 
puted (two significant digits) even for a rather coarse mesh. 



5.2.3.2 Experiment set up based on existing materials 



The method proposed in this section is adapted to z-anisotropic materials, such as transparent 
CaCOs d, LiNbOs [El or Ni:YIG ^ and lossy CoPt or CoPd [22]. Let us now consider 



a trapezoidal (cf. Fig. |5.5| ) anisotropic grating made of aragonite (CaCOs) deposited on an 
isotropic substratum (Si02, £si02 = 2.25). Along the anisotropic crystal axis, its dielectric 
tensor can be written as follows lfT9l : 



=CaC03 




and 



=CaC03 



Mo 
= I Alo 
Mo 



(5.45) 



600 nm 



300 nm 




500 nm 



Fig. 5.5: Diffractive elemen t pattern. This element is made of aragonite for which the dielectric 
tensor is given by Eq. {5.46) and is deposited on a silica substrate with a period d = 600 ^m. 



Now let's assume that the natural axis of our aragonite grating are rotated by 45° around 
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the grating infinite dimension. The dielectric tensor becomes: 



,45° 
=CaC03 




0.251 
2.592 
2.829 



(5.46) 



We shall here remind that our method remains strictly the same whatever the diffractive element 
geometry is. The 2D computational domain is bounded along the };-axis by the PMLs and 
along the x since we consider only one pseudo period. We propose to calculate the diffractive 
efficiencies at Aq = 633 nm in both polarization cases TE and TM, and for different incoming 
incidences (0°, 20° and 40°). Since both ji and £ are Hermitian, the whole incident energy is 
diffracted and the sum of theses efficiencies ought to be equal to the incident energy, which will 
stand for validation of our numerical calculation. 

Finally, the resulting bounded domain is meshed with a maximum mesh element side size 
of Ao/10 v^. Efficiencies are still post-processed in accordance with the calculation presented 
section |5]Z2]6l 



TM 


T-2 


T-x 


7b 


Ti 


R-i 


Rq 


Ri 


total 


0° 




0.203133 


0.585235 


0.203138 




0.008473 




0.999978 


20° 




0.399719 


0.575625 


0.004643 


0.004412 


0.015630 




1.000029 


40° 


0.025047 


0.420714 


0.493491 




0.002541 


0.058238 




1.000031 


TE 


T-2 


T-x 




Ti 


R-l 


Rq 


Rx 


total 


0° 




0.322510 


0.538165 


Q.UAlll 




0.014683 




1.000080 


20° 




0.538727 


0.444403 


0.000369 


0.005372 


0.011180 




1.000051 


40° 


0.012058 


0.434191 


0.541090 




0.005032 


0.007686 




1.000057 



Tab . 5.2 : Transmitted and reflected efficiencies ofpropagative orders deduced from fleld maps shown 



Fig. 5.6 



At normal incidence, the h field in the TE case (cf. Fig. |5.6b| ) is non symmetric whereas 



the e field in the TM ca se is (cf. Fig. |5.6a| ). This is illustrated by the obvious non-symmetry of 
r5f and (cf. Tablels^ 0.322510 versus 0.124722!), whereas T™ = T™ = 0.20313. 
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TM 



TE 



-1.5 -0.5 0.5 1.7 -2.5 -1.5 -0.5 0.5 



1.5 2.5 
^^^^ 





'A 
4 



(a) inV /m at Go = 0° 



(b) %e{h} inA/mat9o = 0° 





(c) inV/m at Oq = 20° 



(d) 'Sie{h} in A/m at do = 20° 




(e) 3ie{e} in V /m at Oq = 40° (f) 3ie{h} in A/m at Oq = 40° 

Fig. 5.6: Real part of the total calculated field depending on Oq and the polarization case 



5.2.3.3 A non trivial geometry 



Since the beginning of this chapter, we have laid great stress upon the independence of the 
method towards the geometry of the pattern. But we have considered so far diffractive objects 
of simple trapezoidal section. Let us tackle a way more challenging case and see what this 
approach is made of. 

We can obtain an quite winding shape by extracting the contrast contour of an arbitrary 
image (see Fig. |5.7^ - |5.7} 3). The contour is approximated by a set of splines, and the resulting 
domain is finely meshed (Fig. 5.7:). Finally, as shown in Fig. 5.7 3, the formed pattern (/^^/Aq = 
1 .68), breathing in free space (^substrat = 1). is supposed to be periodically repeated J/Ao = 1 .26 
on a plane ground of glass (£si02 = 2.25). The element is considered to be "made of a lossy 
material of high optical index (^marsu = 40 + 0.1 /)• The response of this system to a incident 
^--polarized plane wave at oblic incidence (0o = —30°) is finally calculated. The real part of the 



quasi-periodic total field is represented in Fig. |5.7| A for three periods. 

Indeed, we do not have any tabulated data available to check our results. But what we 
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do have is a pretty reliable consistency check through the computation of the energy balance 
described by Eqs. ( |5.42| ) and ( |5.43| ). As shown in Fig. |5.7^ , we obtain at least 7 significative 
digits on the energetic values. The total balance of 1.00000019 is computed taking into account 
(i) values of the total field inside the diffractive elements, (ii) values of the diffracted field at 
altitudes spanning the whole (modeled) superstrate, (iii) values of the total field at altitudes 
spanning the entire (modeled) substrate. Finally, (iv) the calculated field ^2 ^Iso nicely decays 
exponentially inside both PML. These four points allow us to check a posteriori the validity of 
the field everywhere in the computation cell. 





(e) ENERGY BALANCE 



Ro 


0.13436549 




0.10229352 


T-2 


0.00802646 


U 


0.00870891 


To 


0.00666213 




0.00359647 


Q 


0.73634721 


TOTAL 


1.00000019 



Fig. 5.7: (a) Initial contrasted image, (b) Proposed set up. (c) Sample mesh, (d) 3ie{Ez} in V/m. 
(e) Energy balance of the problem. 
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5.2.4 Dealing with Wood anomalies using Adaptative PML 



As we have noticed at the end of Section 5.2.2.4, PMLs based on "traditional coordinate stretch 



ing" are inefficient for periodic problems when dealing with grazing angles of diffracted orders, 
i.e, when the frequency is near a Wood's anomaly ([[23l|24|), leading to spurious reflexions and 
thus numerical pollution of the results. An important question in designing absorbing layers 
is thus the choice of their parameters: The PML thickness and the absorption coefficient. To 
this aim, adaptative formulations have already been set up, most of them employing a posteri- 
ori error estimate t25l[T8l[26l. In this section, we propose Adaptative PMLs (APMLs) with a 
suitable coordinate stretching, depending both on incidence and grating parameters, capable of 
efficiently absorbing propagating waves with nearly grazing angles. This section is dedicated 
to the mathematical formulation used to determine PML parameters adapted to any diffraction 
orders. We provide at the end a numerical example of a dielectric slit grating showing the 
relevance of our approach in comparison with classical PMLs. 

5.2.4.1 Skin depth of the PML 



£^{x,y),jl^{x,y) 




£ {x,y),fL {x,y) 



top PML 



truncated 
superstrate 



groove region 



truncated 
substrate 



bottom PML 



Fig. 5.8: The basic cell used for the FEM computation of the diffracted field 



As explained in Section 



5.2.2.6 



the diffracted field uf^ can be expanded as a Rayleigh 
expansion, i,e, into an infinite sum of propagating and evanescent plane waves called diffraction 
orders. As detailed at the end of Section |5".2.2.4[ we are now in position to rewrite easily the 
expression of, say, a transmitted diffraction order into the substrate. Similar considerations 
also apply to the reflected orders in the top PML. Combining Eq. ( |5.32| ) and ( |5.40[ ) lead to the 
expression u~^{yc) of a transmitted propagative order inside the PML: 



Xyc)^u^ {y{yc)) 



-il5-[yH-{yc-y) 
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The non oscillating part of this function is given by: 

U- (y) = tn exp ((/3^ C"'" + iS^'" C'^yc) , 

where j8~ = Pn~ + i^n~^ For a propagating order we have fin~ > and I5n~ = 0, while for 
an evanescent order Pn~ = and Pn~ > 0. It is thus sufficient to take > and > 
to ensure the exponential decay to zero of the field inside the PML if it was of infinite extent. 
But, of course, for practical purposes, the thickness of the PML is finite and has to be suitably 
chosen. Two pitfalls must be avoided: 

1. The PML thickness is chosen too small compared to the skin depth. As a consequence, 
the electromagnetic wave cannot be considered as vanishing: An incident electromagnetic 
"sees the bottom of the PML". In other words, this PML of finite thickness is no longer 
reflection-less. 

2. The PML thickness is chosen much larger than the skin depth. In that case, a significant 
part of the PML is not useful, which gives rise to the resolution of linear systems of 
unnecessarily large dimensions. 

Then remains to derive the skin depth, /~, associated with the propagating order n. This charac- 
teristic length is defined as the depth below the PML at which the field falls to l/e of its value 
near the surface: 

Finally, we find /" = (j8^ "C" + /3;^' "C")"^ and we define / as the largest value among the 

In 



I = maxL . 



The height of the bottom PML region is set toh =10/" 



5.2.4.2 Weakness of the classical PML for grazing diffracted angles 

Let us consider the (bottom) PML adapted to the substrate. Similar conclusions will hold for the 
top PML. The efficiency of the classical PML fails for grazing diffracted angles, in other words 
when a given order appears/vanishes: this is the so-called Wood's anomaly, well known in the 
grating theory. In mathematical terms, there exists no such that j8~ 0. The skin depth of the 
PML then becomes very large. To compensate this, it is tempting to increase the value of 
but it would lead to spurious numerical reflections due to an overdamping. For a fixed value 
of h~, if is too weak, the absorption in the PMLs is insufficient and the wave is reflected 



on the outward boundary of the PML. To illustrate these typical behaviors (cf. Fig. |5.9| ), we 
compute the field diffracted by a grating with a rectangular cross section of height = 1.5 [im 
and width = 3 |Lim with £^ = 1 1.7, deposited on a substrate with permittivity £~ = 2.25. The 
structure is illuminated by a p-polarized plane wave of wavelength Aq = lOjiim and of angle 
of incidence Gq = W in the air = 1). AH materials are non magnetic (/i^ = 1) and the 
periodicity of the grating is d = 4|Lim. We set h~ = 10/q" and = 1. 
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5.2.4.3 Construction of an adaptative PML 

To overcome the problems pointed out in the previous section, we propose a coordinate stretch- 
ing that rigorously treats the problem of Wood's anomalies. The wavelengths "seen" by the 
system are very different depending on the order at stake: 

• if the diffracted angle Gn is zero, the apparent wavelength Ao / cos Gn is simply the incident 
wavelength, 

• if the diffracted angle is near ±n/2 (grazing angle), the apparent wavelength Aq/ cos Gn 
is very large. 

Thus if a classical PML is adapted to one diffracted order, it will not be for another, and vice 
versa. The idea behind the APML is to deal with each and every order when progressing in the 
absorbing medium. 

Once again the development will be conducted only for the PML adapted to the substrate. 
We consider a real- valued coordinate mapping yd{y), the final complex- valued mapping is then 
ydy) = C~yd{y)^ with the complex constant with > and > 0, accounting for the 
damping of the PML medium. 

We begin with transforming the equation j8^^ = Z:^^ — a^^, so that the function with 
integer argument ^ ^ jS ~ becomes a function with real argument continuously interpolated 
between the imposed integer values. Indeed, the geometric transformations associated to the 
PML has to be continuous and differentiable in order to compute its Jacobian. To that extent, 
we choose the parametrization: 

«(}'d) = ao + — ^, (5.47) 
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so that the apphcation j8~ defined by j8~(_yj)^ = ^qE" — ociy^)^ is continuous. Thus, the prop- 
agation constant of the n^^ transmitted order is given by j8~ = j8~(^Ao). The key idea is to 
combine the complex stretching with a real non uniform contraction (given by the continuous 



function y{yd), Eq. (5.49)). This contraction is chosen in such a way that for each order n there 
is a depth y^ such that, around this depth, the apparent wavelength corresponding to the order 
in play is contracted to a value close to Aq. At that point of the PML, this order is perfectly 
absorbed thanks to the complex stretch. We thus eliminate first the orders with quasi normal 
diffracted angles at lowest depths up to grazing orders (near Wood's anomalies) which are ab- 
sorbed at greater depths. In mathematical words, the translation of previous considerations on 
the real contraction can be expressed as: 

txp[-iP~{yd)y{yd)] =txp{-ikoyd) (5.48) 
The contraction y{yd) is thus given by: 

P [yd) ^y£- -{smOo+yd/dy 



The function y{yd) has two poles, denoted yd±— d{±V £~ — sin 0o)- When y^d±~ =t^^ with 
n G N"*", P~(yd ±) — j8~ (±^Ao) = jS^ = 0, i.e. we are on a Wood's anomaly associated with the 
appearance/disappearance of the ±n^^ transmitted order. We now search for the nearest point to 
^ associated with a Wood's anomaly, denoting: 

nt/ D_= mmJy*_+n*Ao| 

In a second step, we look for the point y^ = n^?io such that: 

fz^/ D= min {D^,D-). (5.50) 

To avoid the singular behaviour at _y j = ±' we continue the graph of the function yd{y) by 
a straight line tangent at y^^, which equation is to{yd) = + where s{yd) = 

■§^{yd) is the so-called stretching coefficient. The final change of coordinate is then given by : 

{y{yd) ioYyd<y^^ 
(5.51) 
to{yd) for yd > yd- 



Figure 5.10 shows an example of this coordinate mapping. Eventually, the complex stretch Sy 



used in Eq. ( |5.29| ) is given by: 

Sy{yd) = C-^{yd)- (5.52) 

oyd 

Equipped with this mathematical formulation, we can tailor a layer that is doubly perfectly 
matched: 

• to a given medium, which is the aim of the PML technique, through Eq. ( |5.27[ ), 

• to all diffraction orders, through the stretching coefficient Sy, which depends on the char- 
acteristics of the incident wave and on opto-geometric parameters of the grating. 
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Fig. 5.10: Example of a coordinate mapping y{yd) used for the APML (black solid line). The graph 
ofyd{y) (blue solid line) is continued by a straight line to{yd) tangent aty^ ( red dashed line) to avoid 
the singular behaviour at y^ = y^. 



5.2.4.4 Numerical example 



We now apply the method described in the preceding parts to design an adapted bottom PML for 



the same example as in part |5. 2.4.2} The parameters are the same, and we choose the wavelength 
of the incident plane wave close to the Wood's anomaly related to the +1 transmitted order 
(Ao = 0.999}^^ ^). Moreover, we set the length of the PML h~ = l.ly^ ^ and choose absorption 
coefficients ^+ = ^~ = 1 + /. For both cases (PML and APML), parameters are alike, the only 
difference being the complex stretch Sy. 



The field maps of the norm of H^, and Ey are plotted in logarithmic scale on Fig. 5.11 



for the case of a classical PML and our APML. We can observe that the field that is effec- 



tively computed is clearly damped in the bottom APML (leftmost on Fig. 5.11 ^b)) whereas it 



is not in the standard case (leftmost on Fig. |5.11| [a)), causing spurious reflections on the outer 
boundary. The fields E^ and Ey are deduced from thanks to Maxwell's equations. The high 



values of Ey at the tip of the APML (rightmost on Fig. 5.11 ^b)) are due to very high values of 



the optical equivalent properties of the APML medium (due to high values of Sy), which does 
not affect the accuracy of the computed field within the domain of interest. 
Another feature of our approach is that it efficiently absorbs the grazing diffraction order, as 
illustrated on Fig. |5.12[ the +1 transmitted order does not decrease in the standard PML (blue 
solid line), and reaches a high value at _y = —hr, whereas the same order tends to zero as 
y ^ —Ir in the case of the adapted PML (blue dashed line). 

To further validate the accuracy of the method, we compare the diffraction efficiencies com- 
puted by our FEM formulation with PML and APML to those obtained by another method. 
We choose the Rigorous Coupled Wave Analysis (RCWA), also known as the Fourier Modal 
Method (FMM, |[27ll ). For the chosen parameters, only the 0^^ order is propagative in reflexion 
and the orders —1,0 and +1 are non evanescent in transmission. We can also check the energy 
balance 5 = 7?o + T^-i + Tq + T+i since there is no lossy medium in our example. Results are 



reported in Table 5.3 , and show a good agreement of the FEM with APML with the results from 
RCWA. On the contrary, if classical PML are used, the diffraction efficiencies are less accurate 
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-2-1 -1 1 2 3 -1 1 2 3 - 3-2-1 - 4-2 2 2 4 




-d/2 d/2 -d/2 t//2 -d/2 C//2 -^//2 d/2 -d/2 i//2 -i//2 d/2 

JC JC JC X JC JC 



(a) Classical PML 



(b) Adapted PML 



Fig. 5.11: Field maps of the logarithm of th e norm o fH^, Ej^ and Eyfor the dielectric slit grating at 
Ao = 0.999jj ^ (same parameters as in part 5.2.4.2). (a): classical PML with inefficient damping of 
Hz in the bottom PML. (b): APML where the field is correctly damped in the bottom PML. For 
both cases the thickness of the PML is h~ = 1.1 ^. 




Fig. 5.12: Modulus of the Un for the three propagating orders with adapted (dashed lines) and 
classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except for the 
grazing one (n= I) as expected. This drawback is bypassed when using the adaptative PML. 



compared to those computed with RCWA. Checking the energy balance leads the same conclu- 
sions: the numerical result is perturbed by the reflection of the waves at the end of the PML if 
it is not adapted to the situation of nearly grazing diffracted orders. 
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Ro 


T-i 


To 


7^1 


B 


RCWA 
FEM + APML 
FEM + PML 


0.1570 
0.1561 
0.1904 


0.3966 
0.3959 
0.4118 


0.1783 
0.1776 
0.1927 


0.2680 
0.2703 
0.2481 


0.9999 
0.9999 
1.0430 



Tab. 5.3: Diffraction efficiencies Rq, T-\, Tq and T+i of the four propagating orders, and energy 
balance B = Rq -\- T-i -\- Tq -\- T^i, computed by three methods: RCWA (line 1), FEM formulation 
with APML (line 2), FEM formulation with classical PML (line 3). 




Fig. 5.13: Mean value of the norm of along the outer boundary of the bottom PML y = 
{\Hz{—h~)\)x,for Ao approaching the Wood's anomaly y'^ ^ by inferior values (Aq = (1 — I0~^)y'^ ^, 
red squares) and by superior value fAo = (1 + I0~^)y^ ^, blue circles) as a function ofn. 



Eventually, to illustrate the behavior of the adaptative PML when the incident wave- 
length gets closer to a given Wood's anomaly, we computed the mean value of the norm of 
Hz along the outer boundary of the bottom PML 7= {\Hz{-h~)\)x, when Aq = (1 + 10"'')}^^ ^ 
and Ao = (1 — I0~^)y^ ^, for ^ = 1,2, ...10. The results are shown in Fig. 5.13 As the wave- 
length gets closer to ^, 7 first increases but for ^ > 3, it decreases exponentially. However, 
in all cases, the value of 7 remains small enough to ensure the efficiency of the PMLs. 



5.2.5 Concluding remarks 

A novel FEM formulation was adapted to the analysis of z-anisotropic gratings relying on a 
rigorous treatment of the plane wave sources problem through an equivalent radiation problem 
with localized sources. The developed approach presents the advantage of being very general 
in the sense that it is applicable to every conceivable grating geometry. 

Numerical experiments based on existing materials at normal and oblique incidences in 
both TE and TM cases showed the efficiency and the accuracy of our method. We demonstrated 
we could generate strongly imbalanced symmetric propagative orders in the TE polarization 
case and at normal incidence with an aragonite grating on a silica substratum. 
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We also introduced the adaptative PML for grazing incidences configurations. It based on 
a complex- valued coordinate stretching that deals with grazing diffracted orders, yielding an ef- 
ficient absorption of the field inside the PML. We provided an example in the TM polarization 
case (but similar results hold for the TE case), illustrating the efficiency of our method. The 
value of the magnetic field on the outward boundary of the PML remains small enough to con- 
sider there is no spurious reflection. The formulation is used with the FEM but can be applied to 
others numerical methods. Moreover, the generalization to the vectorial three-dimensional case 
is straightforward: the recipes given in this last section do work irrespective of the dimension 
and whether the problem is vectorial. 

In the next section, the scalar formulation adapted to mono-dimensional gratings is ex- 
tended to the the most general case of bi-dimensional grating embedded in an arbitrary multi- 
layered dielectric stack with arbitrary incidence. 
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5.3 Diffraction by arbitrary crossed-gratings : a vector Finite Element formulation 
5.3.1 Introduction 



In this section, we extend the method detailed in Sec. |5.2| to the most general case of vector 
diffraction by an arbitrary crossed gratings. The main advantage of the Finite Element Method 
lies in its native ability to handle unstructured meshes, resulting in a build-in accurate discretiza- 
tion of oblique edges. Consequently, our approach remains independent of the shape of the 
diffr active element, whereas other methods require heavy adjustments depending on whether 
the geometry of the groove region presents oblique edges (e.g. RCWA [f28ll . FDTD. . . ). In this 
section, for the sake of clarity, we recall again the rigorous procedure allowing to deal with the 
issue of the plane wave sources through an equivalence of the diffraction problem with a radi- 
ation one whose sources are localized inside the diffractive element itself, as already proposed 
in Sec.|5;2][l9ll3l. 

This approach combined with the use of second order edge elements allowed us to re- 
trieve with a good accuracy the few numerical academic examples found in the literature. Fur- 
thermore, we provide a new reference case combining major difficulties such as a non trivial 
toroidal geometry together with strong losses and a high permittivity contrast. Finally, we dis- 
cuss computation time and convergence as a function of the mesh refinement as well as the 
choice of the direct solver. 



5.3.2 Theoretical developments 

5.3.2.1 Set up of the problem and notations 

We denote by x, y and z the unit vectors of the axes of an orthogonal coordinate system Oxyz. 
We only deal with time-harmonic fields; consequently, electric and magnetic fields are repre- 
sented by the complex vector fields E and H, with a time dependance in txp{—icot). Note that 
incident light is now propagating along the z-axis, whereas j-axis was used in the 2D case. 

Besides, in this section, for the sake of simplicity, the materials are assumed to be isotropic 
and therefore are optically characterized by their relative permittivity £ and relative permeability 
/i (note that the inverse of relative permeabilities are denoted here v). It is of importance to 
note that lossy materials can be studied, the relative permittivity and relative permeability being 
represented by complex valued functions. The crossed-gratings we are dealing with can be split 



into the following regions as suggested in Fig. 5.14 



The superstrate (z > zo) is supposed to be homogeneous, isotropic and lossless, and there- 
fore characterized by its relative permittivity and its relative permeability /i+ (= 1 / v+) 
and we denote ko y^£+/i+, where ko := (o/c, 

The multilayered stack (zn <z<zo) is made of N layers which are supposed to be homo- 
geneous and isotropic, and therefore characterized by their relative permittivity their 
relative permeability /i^(= 1/ v") and their thickness We denote kn := ko for 
n integer between 1 and N. 

The groove region {zg <z< which is embedded in the layer indexed g {e^ of 

the previously described domain, is heterogeneous. Moreover the method does work ir- 
respective of whether the diffractive elements are homogeneous: The permittivity and 
permeability can vary continuously (gradient index gratings) or discontinuously (step 
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index gratings). This region is thus characterized by the scalar fields £^ (x,_y,z) and 
/i^ (x,_y,z)(= 1/v^ {x^y^z}). The groove periodicity along thex-axis, respectively (resp.) 
_y-axis, is denoted dx, resp. dy, in the sequel. 

• The substrate (z < zn) is supposed to be homogeneous and isotropic and therefore char- 
acterized by its relative permittivity £~ and its relative permeability /i~ (= 1/ v~) and we 
denote k~ := ko y^£~/i~, 

Let us emphasize the fact that the method principles remain unchanged in the case of several 
diffractive patterns made of distinct geometry and/or material. 




Fig. 5.14: Scheme and notations of the studied bi-gratings. 

The incident field on this structure is denoted: 

E^^^=A^ exp(/k+.r) 



with 



and 









— sin0 


ocoscpo 


k+ = 


/3o 


= yt+ 


— sin0 






. To . 




— cos 6 


'o 



cos ^^fQ COS 00 COS (po — sin y/b sin cpo 
cos i/Zq cos Qq sin (pQ + sin i/Iq cos (pQ 
— cos \\fQ sin 00 



(5.53) 



(5.54) 
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where (po ^ [0,2;r], 0o ^ [0,7r/2] and i/Aq G [0,;r] (polarization angle). 

We recall here the diffraction problem: finding the solution of Maxwell equations in 
harmonic regime i.e. the unique solution (E,H) of: 

curl E = / 0) ^0 M H (5.56a) 
curlH = -/a)£o^E (5.56b) 

such that the diffracted field satisfies the so-called Outgoing Waves Condition (OWC lISTIl ) and 
where E and H are quasi-bi-periodic functions with respect to x and y coordinates. 



One can choose to calculate arbitrarily E, since H can be deduced from Eq. ( |5.56a| ). The 
diffraction problem amounts to looking for the unique solution E of the so-called vectorial 
Helmholtz propagation equation, deduced from Eqs. ( |5.56a|5.56b| ): 

J^^ ^ := -curl(vcurlE)+/:^£E = (5.57) 

such that the diffracted field satisfies an OWC and where E is a quasi-bi-periodic function with 
respect to x and y coordinates. 



5.3.2.2 From a diffraction problem to a radiative one with localized sources 



According to Fig. |5.14[ the scalar relative permittivity £ and inverse permeability v fields asso- 
ciated to the studied diffractive structure can be written using complex- valued functions defined 
by part and taking into account the notations adopted in Sec. |5.3.2T : 



v{x,y,z) :-- 



v+ 


for 




z>zo 






V" 


for 


Zn- 


-l>Z>Zn 


with 


1 < « < g 


vs'{x,y,z) 


for 




-\>Z>Zg 








for 


Zn- 


-\>Z>Zn 


with 


g <n<N 


v~ 


for 




Z<Zn 







(5.58) 



with V = {£,v} , zo = and = - ^4=1 for I <n< N. 

It is now convenient to introduce two functions defined by part £1 and Vi corresponding to the 
associated multilayered case (i.e. the same stack without any diffractive element) constant over 
Ox and Oy: 

{D+ for z > 
v"" for Zn-i>z>Zn with l<n<N (5.59) 
v~ for z<ZN 

with V = {£, v}. 

We denote by Eq the restriction of E^^^ to the superstrate region: 



E, 



'0 • = 



gmc 





for z > zo 
for z < zo 



(5.60) 



We are now in a position to define more explicitly the vector diffraction problem that we are 
dealing with in this section. It amounts to looking for the unique vector field E solution of: 



^g^(E) = such that E^ : = E - Eq satisfies an OWC. 



(5.61) 



In order to reduce this diffraction problem to a radiation one, an intermediary vector field de- 
noted El is necessary and is defined as the unique solution of: 



(El ) = such that Ef := Ei - Eq satisfies an OWC. 



(5.62) 
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The vector field Ei corresponds to an ancillary problem associated to the general vectorial 
case of a multilayered stack which can be calculated independently. This general calculation is 
seldom treated in the literature, we present a development in Appendix. Thus Ei is from now 
on considered as a known vector field. It is now apropos to introduce the unknown vector field 
E2, simply defined as the difference between E and Ei, which can finally be calculated thanks 
to the FEM and: 

Ef := E - El = E^ - Ef . (5.63) 

It is of importance to note that the presence of the superscript d is not fortuitous: As a difference 
between two diffracted fields (Eq. (5.63 ), E2 satisfies an OWC which is of prime importance in 



our formulation. By taking into account these new definitions, Eq. ( |5.61| ) can be written: 

^,,v(Ef) = -^.,v(Ei), (5.64) 

where the right-hand member is a vector field which can be interpreted as a known vectorial 
source term —S^\{x^y^z) whose support is localized inside the dijfractive element itself. To 
prove it, let us introduce the null term defined in Eq. ( |5.62| ) and make the use of the linearity of 
which leads to: 

^1 :=^,,v(Ei) =^^,v(Ei)-^^i,vi(Ei) =^^-^i,v-vi(Ei) . (5.65) 



5.3.2.3 Quasi-periodicity and weak formulation 

The weak form is obtained by multiplying scalarly Eq. ( 5.61| ) by weighted vectors E^ chosen 
among the ensemble of quasi-bi-periodic vector fields of L^(curl) (denoted (curl, [dx, dy),\i)) 
in Q. '. 

^ev(E,EO= / -curl(vcurlE).E^ + yt^£E.E^d^2 (5.66) 
' Jo. 

Integrating by part Eq. ( |5.66| ) and making the use of the Green-Ostrogradsky theorem lead to: 
^£v(E,EO= / -vcurlE-curlE^ + yt^eE-E^d^l- / (n x (v curlE)) -E^dS (5.67) 

where n refers to the exterior unit vector normal to the surface dQ. enclosing Q.. 

The first term of this sum concerns the volume behavior of the unknown vector field 
whereas the right-hand term can be used to set boundary conditions (Dirichlet, Neumann or so 
called quasi-periodic Bloch-Floquet conditions). 

The solution E2 of the weak form associated to th e diffraction problem, expressed in its 
previously defined equivalent radiative form at Eq. (5.64), is the element of (curl, {dx, dy),\i) 
such that: 



VE^ GL2(curl,4,rf>.,k),^,,v(Ef,E0 = -^^-^i,v-vi(Ei,EO . (5.68) 

In order to rigorously truncate the computation a set of Bloch boundary conditions are 
imposed on the pair of planes defined by {y — —dyjl^y — dy/2) and (x = —dx/2.^x — dx/2). 
One can refer to [1 1 J for a detailed implementation of Bloch conditions adapted to the FEM. A 
set of Perfectly Matched Layers are used in order to truncate the substrate and the superstrate 
along z axis (see [[32l for practical implementation of PML adapted to the FEM). Since the 
proposed unknown E2 is quasi-bi-periodic and satisfies an OWC, this set of boundary conditions 
is perfectly reasonable: E2 is radiated from the diffractive element towards the infinite regions 
of the problem and decays exponentially inside the PMLs along z axis. The total field associated 
to the diffraction problem E is deduced at once from Eq. ( |5.63| ). 
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In the vectorial case, edge elements (or Whitney forms) make a much more relevant choice ff 33]| 
than nodal elements. Note that a lot of work (see for instance [l34l|) has been done on higher 
order edge elements since their introduction by Bossavit ll35]| . These elements are suitable to the 
representation of vector fields such as by letting their normal component be discontinuous 
and imposing the continuity of their tangential components. Instead of linking the Degrees Of 
Freedom (DOF) of the final algebraic system to the nodes of the mesh, the DOF associated to 
edges (resp. faces) elements are the circulations (resp. flux) of the unknown vector field along 
(resp. across) its edges {rt^^. faces). 

Let us consider the computation cell together with its exterior boundary dO.. This 
volume is sampled in a finite number of tetrahedron according to the following rules: Two 
distinct tetrahedrons have to either share a node, an edge or a face or have no contact. Let us 
denote by the set of tetrahedrons, ^ the set of faces, the set of edges and J/^ the set of 
nodes. In the sequel, one will refers to the node n — {/}, the edge e — {/, j}, the face / = {/, k} 
and the tetrahedron t — {/, j, /}. 




Fig. 5.15: Degrees of freedom of a second order tetrahedral element. 



Twelve DOF (two for each of the six edges of a tetrahedron) are classically derived from 
line integral of weighted projection of the field E2 on each oriented edge e — {ij}'. 

f^Ei-tijXidi 

, (5.69) 

= J^Ei-tjiXjdl 

where t/j is the unit vector and A/, the bary centric coordinate of node /, is the chosen weight 
function. 

According to Yioultsis et al. Il36l . a judicious choice for the remaining DOF is to make 
the use of a tangential projection of the 1-form E2 on the face / = {/, j, fc}. 



^ikj = J X n,::^) • gradA;td5 



(5.70) 
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The expressions for the shape functions, or basis vectors, of the second order 1-form Whitney 
element are given by: 



: (8 A^^ - 4 A/) grad A^- + (- 8 A/ Xj + 2 A^-) grad A/ 
: 1 6 A/ A j grad Ay^ — 8 Xj grad A/ — 8 Ay^ A/ grad Ay 



(5.71) 



This choice of shape function ensures OTll the following fundamental property: every degree of 
freedom associated with a shape function should be zero for any other shape function. Finally, 



an approximation of the unknown E2 projected on the shape functions of the mesh m (E2 
be derived: 

-^d,m 



) can 



(5.72) 



Weight functions E^ (c.f. Eq. (5.68) are chosen in the same space than the unknown E2, 
L^(curl, {dx,dy),k). According to the Galerkin formulation, this choice is made so that their 
restriction to one bi-period belo ngs t o the se t of s hape functions mentioned above. Inserting 
the decomposition of E2 of Eq. (5.72) in Eq. (5.68 ) leads to the final algebraic system which is 
solved, in the following numerical examples, thanks to direct solvers. 



5.3.3 Energetic considerations: Diffraction efficiencies and losses 

Contrarily to modal methods based on the determination of Rayleigh coefficients, the rough 
results of the FEM are three complex components of the vector field interpolated over the 
mesh of the computation cell. Diffraction efficiencies are deduced from this field maps as 
follows. 

As a difference between two quasi-periodic vector fields (see Eq. ( 5.61 )), E^ is quasi-bi- 
periodic and its components can be expanded as a double Rayleigh sum: 



(f7,m)GZ2 



with = ao + jS^ = jSo + and 



(5.73) 



d^x 



iz) 



I rd,/2 rdy/2 

j-d^/2 j-dy/2 



d 



(5.74) 



By inserting the decomposition of Eq. (5.73), which is satisfied by everywhere but in the 
groove region, into the Helmholtz propagation equation, one can express Rayleigh coefficients 
in the substrate and the superstrate as follows: 



,d,x 



(5.75) 



with = k^^ — — j8^, where jn^m (or —iyn.m) is positive. The quantity ui'^ is the sum of 
a propagative plane wave (which propagates towards decreasing values of z, superscript p) and 
of a counterpropagative one (superscript c). The OWC verified by E^ imposes: 



V(fz,m) gZ^ 



x^p 
x,c 



: for z > zo 
for z < za^ 



(5.76) 



Eq. (5.74) allows to evaluate numerically (resp. Cn^) by double trapezoidal integration of 
a slice of the complex component E^ at an altitude Zc fixed in the superstrate (resp. substrate). 
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It is well known that the mere trapezoidal integration method is very efficient for smooth and 
periodic functions (integration on one period). The same holds for Ey and E^ components as 

well as their coefficients en}m^^ and enjn'^^- 

The dimensionless expression of the efficiency of each reflected and transmitted {n,m) 
order [381 is deduced from Eqs. ( |5.75|5.76| ): 

1 ?n,m c 



" ^0 (5.77) 



5 



Tfi^m — A2 ' ^n,m{Zc) ' ^n,m{Zc) for Zc < Zn 

with e^^'^^ - /'^^'^^ X + /'^^'^^ V + /'^^'^^ z 

Furthermore, normalized losses Q can be obtained through the computation of the follow- 
ing ratio: 

Q = — . (5.78) 

-9?4EoxHo}-nd5 



Is 2 

The numerator in Eq. ( |5.78| ) clarifies losses in watts by bi-period of the considered crossed- 
grating and are computed by integrating the Joule effect losses density over the volume V of 
the lossy element. The denominator normalizes these losses to the incident power, i.e. the time- 
averaged incident Poynting vector flux across one bi-period (a rectangular surface S of area 
dxdy in the superstrate parallel to Oxy, whose normal oriented along decreasing values of z is 
denoted n). Since Eq is nothing but the plane wave defined at Eqs. ( |5.54|5.55l ), this last term 



is equal to (A^ ^£o/Mo^x^j)/(2cos(0o))- Volumes and normal to surfaces being explicitly 
defined, normalized losses losses Q are quickly computed once E determined and interpolated 
between mesh nodes. 

Finally, the accuracy and self-consistency of the whole calculation can be evaluated by 
summing the real part of transmitted and reflected efficiencies {n^m) to normalized losses: 

Q+ ^e{Rn,m}+ M7;,m}, 

quantity to be compared to 1 . The sole diffraction orders taken into account in this conservation 
criterium correspond to propagative orders whose efficiencies have a non-null real part. Indeed, 
diffraction efficiencies of evanescent orders, corresponding to pure imaginary values of for 



higher values of {n,m) (see Eq. (5.75)) are also pure imaginary values as it appears clearly in 



Eq. ( |5.77| ). Numerical illustrations of such global energy balances are presented in the next 



section. 
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5.3.4 Accuracy and convergence 

5.3.4.1 Classical crossed gratings 

There are only a few references in the hterature containing numerical examples. For each of 
them, the problem only consists of three regions (superstrate, grooves and substrate) as summed 
up on Figure |5.16[ For the four selected cases, among six found in the literature, published 




Fig. 5.16: Configuration of the studied cases. 

results are compared to ones given by our formulation of the FEM. Moreover, in each case, a 
satisfying global energy balance is detailed. Finally a new validation case combining all the 
difficulties encountered when modeling crossed-gratings is proposed: A non-trivial geometry 
for the diffractive pattern (a torus), made of an arbitrary lossy material leading to a large step 
of index and illuminated by a plane wave with an oblique incidence. Convergence of the FEM 



calculation as well as computation time will be discussed in Sec. 5.3.4.2 
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Checkerboard grating In this example worked out by L. Li ll27ll . the diffractive element 



5.16 



is a recta ngular parallelepiped as shown Fig. |5.17a| and the grating parameter highlighted in 
Fig. 



are the following: (pQ = do = 0°, y/b = 45°, dx = dy = 5 Aq \/2/4, h = Xo,e^ = = 
2.25 and £-=£^ = 1. 




(a) (b) 
Fig. 5.17: Diffractive element with vertical edges (a). ^e{Ex} in V/m (b). 





FMM (TT\ 


FEM 




0.04308 


0.04333 


T-i,o 


0.12860 


0.12845 


T-1,+1 


0.06196 


0.06176 


7b,-i 


0.12860 


0.12838 


7b,o 


0.17486 


0.17577 


^0,+l 


0.12860 


0.12839 




0.06196 


0.06177 


7+1,0 


0.12860 


0.12843 


T+\,+i 


0.04308 


0.04332 


^e{Rn,„} 




0.10040 








TOTAL 




1.00000 



Tab. 5.4: Energy balance ^27\l . 



Our formulation of the FEM shows good agreement with the Fourier Modal Method de- 
veloped by L. Li ( GTll, 1997) since the maximal relative difference between the array of values 
presented in Table |5.4| remains lower than 10"^. Moreover, the sum of the efficiencies of prop- 
agative orders given by the FEM is very close to 1 in spite of the addition of all errors of 
determination upon the efficiencies. 
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Pyramidal crossed-grating In this example firstly worked out by Derrick et al [[39l, the 
diffractive element is a pyramid with rectangular basis as shown Fig. |5.18a| and the grating 
parameters highlighted in Fig. |5.16| are the following: Aq = 1.533, (jpo = 45°, 0o = 30°, i/Aq = 0°, 
4 = 1 .5, J^; = 1, /i = 0.25, = 8^ = 1 and £~ = = 2.25. Results given by the FEM show 




(a) (b) 

Fig. 5.18: Diffractive element with oblique edges (a). ^e{Ey} in V/m (b). 
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7b,-i 
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0.00692 
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0.00303 


0.00370 


0.00294 


0.00294 


0.00299 


7b,o 


0.96219 


0.96316 


0.96472 


0.96448 


0.96447 


7^1,0 


0.00299 


0.00332 


0.00280 


0.00282 


0.00290 


TOTAL 


0.99855 


1.00001 


1.00008 


0.99999 


1.00004 



Tab. 5.5: Comparison with the results given in l(39\\4Ul\4l\\42\l . 



good agreement with ones of the C method [^|43, the Rayleigh method gOl and the RCWA 
ll4T1l . Note that, in this case, some edges of the diffractive element are oblique. 
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Bi-sinusoidal grating In this example worked out by Bruno et al. [|43]| . the surface of the 
grating is bi-sinusoidal (see Fig. |5. 19a] ) and described by the function / defined by: 



f{x.y) 



27ix\ /2;r_y 
cos I I + cos 



(5.79) 



The grating parameters et a/. highlighted in Fig. |5.16| are the following: Aq = 0.83, (po = = 
i/Ao = 0°, dx = dy = l,h = 0.2, 8+ = £^ = 1 and £~ = = 4. Note that in order to define this 





(a) (b) 
Fig. 5.19: Dijfractive element with oblique edges (a). ^e{E^} in V/m (b). 
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TOTAL 




0.99806 



Tab. 5.6: Energy balance l[43\l . 



surface, the bi-sinusoid was first sampled (15x15 points), then converted to a 3D file format. 
This sampling can account for the slight differences with the results obtained using the method 
of variation of boundaries developed by Bruno et at. (1993). 
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Circular apertures in a lossy layer In this example worked out by Schuster et al. O, the 
diffractive element is a circular aperture in a lossy layer as shown Fig. |5.20a| and the grating 
parameter highlighted in Fig. 5.16 are the following: Aq = 500 nm, (jpQ = 0o = 0°, £^ = £^ = 1, 
=0.8125 + 5.2500/ and =2.25. 




(a) (b) 



Fig. 5.20: Lossy diffractive element with vertical edges (a). ^e{Ey} in V/m (b). 
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fT7] 


(445 


FEM 
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0.29110 
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0.44148 


TOTAL 








1.00019 



Tab. 5.7: Comparison with l{45\ \27\ \44\l and energy balance. 



In this lossy case, results obtained with the FEM show good agreement with the ones 
obtained with the FMM [27J, the differential method [4411461 and the RCWA [45 J. Joule losses 
inside the diffractive element can be easily calculated, which allows to provide a global energy 
balance for this configuration. Finally, the convergence of the value 7?o,o as a function of the 
mesh refinement will be examined. 
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Lossy tori grating We finally propose a new test case for crossed-grating numerical methods. 
The major difficulty of this case lies both in the non trivial geometry (see Fig. |5.21a| ) of the 
diffractive object and in the fact that it is made of a material chosen so that losses are optimal 



inside it. The grating parameters highlighted in Fig. |5.16| and Fig. |5.21a| are the following 
Ao = l, (po = V^o = 0°, 4 
£^' = -21+20/and£- = 



J - 0.3, a = 0.1, b = 0.05, R = 0.15, h = 500nm, = £^ = 1, 



7 

2.25. 




(a) (b) 

Fig. 5.21: Torus parameters (a). Coarse mesh of the computational domain (b). 



FEM 3D 


0=0° 


= 40° 


^0,0 


0.36376 
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0.32992 


0.38191 


Q 


0.30639 


0.34476 


TOTAL 


1.00007 


0.99998 



Tab. 5.8: Energy balances at normal and oblique incidence. 



Tab. |5.8| illustrates the independence of our method towards the geometry of the diffractive 
element. £^ is chosen so that the skin depth has the same order of magnitude as b, which max- 
imizes losses. Note that energy balances remain very accurate at normal and oblique incidence, 
in spite of both the non-triviality of the geometry and the strong losses. 
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5.3.4.2 Convergence and computation time 



Convergence as a function of mesh refinement When using modal methods such as the 
RCWA or the differential method, based on the calculation of Rayleigh coefficients, a number 
proportional to A^^^ have to be be determined a priori. Then, the unknown diffracted field is 
expanded as a Fourier serie, injected under this form in Maxwell equations, which annihilates 
X— and J— dependencies. This leads to a system of coupled partial differential equations whose 
coefficients can structured in a matrix formalism. The resulting matrix is sometimes directly 
invertible (RCWA) depending on whether the geometry allows to suppress the z— dependance, 
which makes this method adapted to diffractive elements with vertically (or decomposed in 
staircase functions) shaped edge. In some other cases, one has to make the use of integral 
methods in order to solve the system, as in the pyramidal case for instance, which leads to the so- 
called differential method. The diffracted field map can be deduced from these coefficients. If 
the grating configuration only calls for a few propagative orders and if the field inside the groove 
region is not the main information sought for, these two close methods allow to determine the 
repartition of the incident energy quickly. However, if the field inside the groove region is the 
main piece of information, it is advisable to calculate many Rayleigh coefficients corresponding 
to evanescent waves which increases the computation time as (Nr)^ or even (Nr)^. 

FEM relies on the direct calculation of the vectorial components of the complex field. 
Rayleigh coefficients are determined a posteriori. The parameter limiting the computation time 
is the number of tetrahedral elements along which the computational domain is split up. We 



suppose that it is necessary to calcu 
of the field (Xo/V^e{£}). Figure 



ate a t least two or three points (or mesh nodes) per period 



5.22 shows the convergence of the efficiency 7?o,o (circu- 



lar apertures case, see Fig. 5.20a) as a function of the mesh refinement characterized by the 



parameter Nm- The maximum size of each element is set to Xq/{Nm ^/'^e{£}). 
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Fig. 5.22: Convergence ofRo^o in function ofNm (circular apertures crossed-grating). 
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It is of interest to note that even if Nm < 3 the FEM still gives pertinent diffraction effi- 
ciencies: 7?o,o = 0.2334 for Nm = 1 and /?o,o = 0.2331 for Nm = 2. The Galerkin method (see 
Eq. ( |5.67| )) corresponds to a minimization of the error (between the exact solution and the ap- 
proximation) with respect to a norm that can be physically interpreted in terms of energy-related 
quantities. Therefore, the finite element methods usually provide energy-related quantities that 
are more accurate than the local values of the fields themselves. 



Computation time All the calculations were performed on a server equipped with 8 dual 
core Itaniuml processors and 256Go of RAM. Tetrahedral quadratic edge elements were used 
together with the direct solver PARDISO. Among different direct solvers adapted to sparse 
matrix algebra (UMFPACK, SPOOLES and PARDISO), PARDISO turned out to be the less 



time-consuming one as shown in Tab 5.9 



Solver 


Computation time for 41720 DOF 


Computation time for 205198 DOF 


SPOOLES 


15mn32s 


14h44mn 


UMFPACK 


2mn07s 


lhl2mn 


PARDISO 


57 s 


16mn 



Tab. 5.9: Computation time variations from solver to solver 




Fig. 5.23: Computation time and number of DOF as a function ofNu- 



Figure 5.23 shows the computation time required to perform the whole FEM computa- 
tional process for a system made of a number of DOF indicated on the right-hand ordinate. 
It is of importance to note that for values of Nm lower than 3, the problem can be solved in 
less than a minute on a standard laptop (4Go RAM, 2 x 2 GHz) with 3 significant digits on the 
diffraction efficiencies. This accuracy is more than sufficient in numerous experimental cases. 
Furthermore, as far as integrated values are at stake, relatively coarse meshes {Nm ^1) can be 
used trustfully, authorizing fast geometric, spectral or polarization studies. 

Nowadays, the efficiency of the numerical algorithms for sparse matrix algebra together 
with the available power of computers and the fact that the problem reduces to a basic cell with a 
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size of a small number of wavelengths make the finite element problem very tractable as proved 
here. 

5.4 Concluding remarks 

In this chapter, we demonstrate a general formulation of the FEM allowing to calculate the 
diffraction efficiencies from the electromagnetic field diffracted by arbitrarily shaped gratings 
embedded in a multilayered stack lightened by a plane wave of arbitrary incidence and polar- 
ization angle. It relies on a rigorous treatment of the plane wave sources problem through an 
equivalent radiation problem with localized sources. Bloch conditions and a new dedicated 
PML have been implemented in order to rigorously truncate the computational domain. 

The principles of the method were discussed in detail for mono-dimensional gratings 
in TE/TM polarization cases (2D or scalar case) in a first part, and for the most general bi- 
dimensional or crossed gratings (3D or vector case) in a second part. Note that the very same 
concepts could be applied to the intermediate case of mono-dimensional gratings enlighten by 
an arbitrary incident plane wave (so-called conical case). The reader will find detail about the 
element basis relevant to this case in [UTll . 

The main advantage of this formulation is its complete generality with respect to the 
studied geometries and the material properties, as illustrated with the lossy tori grating non- 
trivial case. Its principle remains independent of both the number of diffractive elements by 
period and number of stack layers. Its flexibility allowed us to retrieve with accuracy the few 
numerical academic examples found in the literature and established with independent methods. 

The remarkable accuracy observed in the case of coarse meshes, makes it a fast tool for 
the design and optimization of diffractive optical components (e.g. reflection and transmission 
filters, polarizers, beam shapers, pulse compression gratings. . . ). The complete independence of 
the presented approach towards both the geometry and the isotropic constituent materials of the 
diffractive elements makes it a handy and powerful tool for the study of metamaterials, finite- 
size photonic crystals, periodic plasmonic structures. . . The method described in this chapter 
has already been successfully applied to various problems, from homogenization theory [47J 
or transformation optics [48] to more applied concerns as the modeling of complex CMOS 
nanophotonic devices ll49l or ultra- thin new generation solar cells lISOll . 



G. Demesy et ah: Finite Element Method 



5.41 



5.A APPENDIX 

This appendix is dedicated to the determination of the vector electric field in a dielectric stack 
enlightened by a plane wave of arbitrary polarization and incidence angle. This calculation, 
abundantly treated in the 2D scalar case, is generally not presented in the literature since, as 
far as isotropic cases are concerned, it is possible to project the general vectorial case on the 
two reference TE and TM cases. However, the presented formulation can be extended to a fully 
anisotropic case for which this TE/TM decoupling is no longer valid and the three components 
of the field have to be calculated as follows. 



Let us consider the ancillary problem mentioned in Sec. |5.3.2.2[ i.e. a dielectric stack 
made of N homogeneous, isotropic, lossy layers characterized by there relative permittivity 
denoted and their thickness ej. This stack is deposited on a homogeneous, isotropic, possibly 
lossy substrate characterized by its relative permittivity denoted £^^^ = The superstrate is 
air and its relative permittivity is denoted = 1. Finally, we denote by Zj the altitude of 
the interface between the and j layers. The restriction of the incident field E^^^ to 
the superstrate region is denoted Eq. The problem amounts to looking for (Ei,Hi) satisfying 
Maxwell equations in harmonic regime (see Eqs. ( |5.56a|5.56b| )). 



Across the interface z = zj 

By projection on the main axis of the vectorial Helmholtz propagation equation (Eq. ( |5.57| )), 
the total electric field inside the layer can be written as the sum of a propagative and a 
counter-propagative plane waves: 



Ei(x,};,z) 



where 



exp{j{aox + poy + rjz)) + 



exp(;(aox + ^oy-r;2)) 
(5.80) 

(5.81) 



What follows consists in writing the continuity of the tangential components of (Ei ,Hi) across 
the interface z = zj, i.e. the continuity of the vector field ^ defined by: 



Ef 
El 



(5.82) 



The continuity of ^ along Oz together with its analytical expression inside the and j + 
layers allows to establish a recurrence relation for the interface z = zj. 
Then, by projection of Eqs. ( |5.56a|5.56bl ) on Ox,Oy and Oz'. 





' El ' 


= -io)e 


El 




.E\. 



(5.83) 
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and 



'1 dz 



Hf 
Hi 



(5.84) 



^-i(^E\ 
_ iaoEl — ipoEf 

Consequently, tangential components of Hi can be expressed in function of tangential compo- 
nents of El! 

^coii jSo 



(Oil -Oq 



iHf 



dEl 

dz 
dEf 

dz 





(5.85) 



By noticing the invariance and linearity of the problem along Ox and Oy, the following notations 
are adopted: 



j,± 



and 



£^'^'2exp(±/7;z) 
if-'' exp(±/7yz) 

Uy 



Thanks to Eq. (5.80) and Eq. (5.84) and letting M = it comes for the layer: 





" 1 


1 





" 




\ 1 


exp(/(aox + /3oy)) 








1 


1 






























TJ-J 



n,- 



Finally, the continuity of ^ at the interface z = Zj leads to: 

^;+i(^i) = n7^in,-4>Xz,-). 

Normal components can be deduced using Eqs. ( |5.83|5. 



(5.86) 



(5.87) 



(5.88) 



(5.89) 



Traveling inside the j+V^ layer 



Using Eq. (5.80), a simple phase shift allows to travel from z = Zj to z = zj+i = Zj — ej+V- 



^;+i(^y+i) = 



exp(-/7y+ie^-+i) 

exp(+/7,-+i6';+i) 








exp(-/7,-+i6'^-+i) 

exp(+/7^-+ie^-+i)_ 



^7+1 



(5.90) 

Thanks to Eq. ( |5.90| ) and Eq. ( |5.89[ ), a recurrence relation can be formulated for the analytical 
expression of Ei in each layer: 



^j+i{zj+x) = Tj+^n-:l,nj^j{zj) 



(5.91) 
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Reflection and transmission coefficients 



The last step consists in the determination of the first term Oq, which is not entirely known, 

imposed by the incident field Eq. 

1 



since the problem definition only specifies Ux'^ and Uy' 

Let us make the use of the OWC hypothesis verified by Ef (see Eq. (15. 621)). This hypothesis 



directly translates the fact that none of the components of Ef can either be traveling down in 



the superstrate or up in the substrate: 



0. 



Therefore, the four unknowns 



Uy' ,Uy^^'^ and Ux^^'^, i-e. transverse components of the vector fields reflected and 



transmitted by the stack, verify the following equation system: 



(5.92) 



This allows to extend the definition of transmission and reflection widely used in the scalar case. 
Fina lly, O a^+i is entirely defined. Making the use of the recurrence relation of Eq. ( |5.91| ) and of 
Eq. (5.80) leads to an analytical expression for Ef in each layer. 
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